
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Empirical Application(Timber Dataset)                                     %
%  Code used to estimate:                                                    %
%  (i) private values quantile regression coefficients;                      %
%  (ii) optimal screening level;                                             %
%  (iii) optimal reservation price;                                          %
%  (iv) Seller's expected  payoff under an optimal reservation price policy; %
%  (v) 95% Confidence Intervals of (i)-(iv) computed via pairwise bootstrap. %             
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all

load('data.mat'); % load timber dataset
L2=size(wb2,1); % number of auctions with N=2
L3=size(wb3,1); % number of auctions with N=3
L=L2+L3;
N=[2,3]; % Number of bidders
alpha=0.10:0.02:0.92; % quantile levels
psi2=ones(L2,1)*(N(1)*(alpha.^(N(1)-1))-(N(1)-1)*(alpha.^N(1))); % Psi function given N=2
psi3=ones(L3,1)*(N(2)*(alpha.^(N(2)-1))-(N(2)-1)*(alpha.^N(2))); % Psi function given N=3
psi = [psi2;psi3];
m=size(alpha,2);

% covariates
x2=[ones(L2,1),apv2,vol2]; 
x3=[ones(L3,1),apv3,vol3]; 
x = [x2;x3];
xbar=quantile(x,0.5); % median auction (it is also considered the quantiles 0.15 and 0.85 of x)
p=size(x,2);

% winning bids
wb=[wb2;wb3];
wbx2=[wb2,x2];   
wbx3=[wb3,x3]; 

% Seller's Private Value
V02=zeros(L2,1); % the other option is V02=apv2
V03=zeros(L3,1); % the other option is V03=apv3
V0=[V02;V03];
V0bar=quantile(V0,0.5);

B=5000; % Bootstrap replications
toler=0.01;  % Parameter of tolerance for precision of the estimation

gammai=ones(p,1); % Initial guess for the parameters

for k=1:m
        [~,~,g]=MM(psi(:,k), p, wb, x, toler, gammai); % estimation of the coefficients
        gk(:,k)=g;
        V_adjbar=xbar*g; % Adjusted private value quantile
        
        %Computation of the the seller expected payoff
        ER1_N2(k)=V_adjbar*(N(1)*(N(1)-1)*(alpha(k)^(N(1)-2))*(1-alpha(k))); % estimate of the integrand in eq (2.9) for each k given N=2
        ER1_N3(k)=V_adjbar*(N(2)*(N(2)-1)*(alpha(k)^(N(2)-2))*(1-alpha(k))); % estimate of the integrand in eq (2.9) for each k given N=3
        ER2_N2(k)=V0bar*alpha(k)^N(1) + V_adjbar*N(1)*(alpha(k)^(N(1)-1))*(1-alpha(k)); % estimate of the first two terms in eq (2.9) for each k given N=2
        ER2_N3(k)=V0bar*alpha(k)^N(2) + V_adjbar*N(2)*(alpha(k)^(N(2)-1))*(1-alpha(k)); % estimate of the first two terms in eq (2.9) for each k given N=3
        
end

gammas=reshape(gk,p*m,1);

% Computation of expected payoff when screening level is in [alpha(1),alpha(m-1)]
for k=1:m-1  
    intmed_N2=trapz(alpha(k:m),ER1_N2(k:m),2); % numerical integration of ER1_N2 from alpha(k) to alpha(m-1) using a trapezoidal rule
    intmed_N3=trapz(alpha(k:m),ER1_N3(k:m),2); % numerical integration of ER1_N3 from alpha(k) to alpha(m-1) using a trapezoidal rule
    ER_N2(:,k)=ER2_N2(k)+intmed_N2; % estimate of the expected payoff given N=2
    ER_N3(:,k)=ER2_N3(k)+intmed_N3; % estimate of the expected payoff given N=3    
end
        
% Expected payoff all alpha in [alpha(1),alpha(m)]
ER_N2all=[ER_N2,ER2_N2(m)]; % estimate given N=2
ER_N3all=[ER_N3,ER2_N3(m)]; % estimate given N=3
        
% Maximization of expected payoff N=2
ER_N2max=max(ER_N2all,[],2); % max of expected payoff
[~,indx]=ismember(ER_N2max,ER_N2all);
qindx_N2=indx; % index in alpha that maximize estimated expected payoff
q_rstar_N2=alpha(qindx_N2); % optimal screening level
psi_q_rstar_N2_N2=N(1)*(q_rstar_N2^(N(1)-1))-(N(1)-1)*(q_rstar_N2^N(1)); % psi function used in the computation of the optimal reservation price
psi_q_rstar_N2_N3=N(2)*(q_rstar_N2^(N(2)-1))-(N(2)-1)*(q_rstar_N2^N(2)); % psi function used in the computation of the optimal reservation price
psi_q_rstar_N2 = [ones(L2,1)*psi_q_rstar_N2_N2;ones(L3,1)*psi_q_rstar_N2_N3];
[~,~,g_q_rstar_N2]=MM(psi_q_rstar_N2, p, wb, x, toler, gammai); % gamma_hat corresponding to the optimal screening level
% Optimal Res Price
V_rstar_N2=xbar*g_q_rstar_N2; % optimal reservation price

% Maximization of expected payoff N=3
ER_N3max=max(ER_N3all,[],2); % max of expected payoff
[~,indx]=ismember(ER_N3max,ER_N3all);
qindx_N3=indx; % index in alpha that maximize estimated expected payoff
q_rstar_N3=alpha(qindx_N3); % optimal screening level
psi_q_rstar_N3_N2=N(1)*(q_rstar_N3^(N(1)-1))-(N(1)-1)*(q_rstar_N3^N(1)); % psi function used in the computation of the optimal reservation price
psi_q_rstar_N3_N3=N(2)*(q_rstar_N3^(N(2)-1))-(N(2)-1)*(q_rstar_N3^N(2)); % psi function used in the computation of the optimal reservation price
psi_q_rstar_N3 = [ones(L2,1)*psi_q_rstar_N3_N2;ones(L3,1)*psi_q_rstar_N3_N3];
[~,~,g_q_rstar_N3]=MM(psi_q_rstar_N3, p, wb, x, toler, gammai); % gamma_hat corresponding to the optimal screening level
% Optimal Res Price
V_rstar_N3=xbar*g_q_rstar_N3; % optimal reservation price

% Computation of the seller expected payoff under optimal screening level
ER2_N2_rstar=V0bar*q_rstar_N2^N(1) + (V_rstar_N2*N(1))*(q_rstar_N2^(N(1)-1))*(1-q_rstar_N2);
ER2_N3_rstar=V0bar*q_rstar_N3^N(2) + (V_rstar_N3*N(2))*(q_rstar_N3^(N(2)-1))*(1-q_rstar_N3);
q2_r=linspace(q_rstar_N2,0.92,m);
q3_r=linspace(q_rstar_N3,0.92,m);
psiN2_q2=N(1)*(q2_r.^(N(1)-1))-(N(1)-1)*(q2_r.^N(1));
psiN3_q2=N(2)*(q2_r.^(N(2)-1))-(N(2)-1)*(q2_r.^N(2));
psiN2_q3=N(1)*(q3_r.^(N(1)-1))-(N(1)-1)*(q3_r.^N(1));
psiN3_q3=N(2)*(q3_r.^(N(2)-1))-(N(2)-1)*(q3_r.^N(2));

for k=1:m
    [~,~,g2_r]=MM([ones(L2,1)*psiN2_q2(k);ones(L3,1)*psiN3_q2(k)], p, wb, x, toler, gammai);
    ER1_N2_rstar(:,k)=xbar*g2_r*(N(1)*(N(1)-1)*(q2_r(k)^(N(1)-2))*(1-q2_r(k)));
    [~,~,g3_r]=MM([ones(L2,1)*psiN2_q3(k);ones(L3,1)*psiN3_q3(k)], p, wb, x, toler, gammai);
    ER1_N3_rstar(:,k)=xbar*g3_r*(N(2)*(N(2)-1)*(q3_r(k)^(N(2)-2))*(1-q3_r(k)));
end
int2_rstar=trapz(q2_r,ER1_N2_rstar,2); % numerical integration of ER1_N2_rstar using a trapezoidal rule
int3_rstar=trapz(q3_r,ER1_N3_rstar,2); % numerical integration of ER1_N2_rstar using a trapezoidal rule

% Expected payoff under optimal reservation price policy 
ER_N2_rstar=ER2_N2_rstar+int2_rstar;
ER_N3_rstar=ER2_N3_rstar+int3_rstar;

% Pairwise Bootstrap to compute confidence intervals
for b=1:B
    
    wbxb2 = wbx2(randsample(1:L2,L2,true),:); % Resampling data in the subsample N=2
    wbxb3 = wbx3(randsample(1:L3,L3,true),:); % Resampling data in the subsample N=3
    wbxb=[wbxb2;wbxb3];
    wb = wbxb(:,1);
    x = wbxb(:,2:p+1);
    xbar=quantile(x,0.5); % median auction (it is also considered the quantiles 0.15 and 0.85 of x)
    
    % Seller's Private Value
    V02=zeros(L2,1); % the other option is V02=apv2
    V03=zeros(L3,1);
    V0=[V02;V03];
    V0bar=quantile(V0,0.5);
    
     for k=1:m
         [~,~,gb]=MM(psi(:,k), p, wb, x, toler, gammai); 
         gkb(:,k)=gb;
         V_adjbar=xbar*gb; % Adjusted private value quantile
         
         %Computation of the the seller expected payoff
         ER1_N2(k)=V_adjbar*(N(1)*(N(1)-1)*(alpha(k)^(N(1)-2))*(1-alpha(k))); % estimate of the integrand in eq (2.9) for each k given N=2
         ER1_N3(k)=V_adjbar*(N(2)*(N(2)-1)*(alpha(k)^(N(2)-2))*(1-alpha(k))); % estimate of the integrand in eq (2.9) for each k given N=3
         ER2_N2(k)=V0bar*alpha(k)^N(1) + V_adjbar*N(1)*(alpha(k)^(N(1)-1))*(1-alpha(k)); % estimate of the first two terms in eq (2.9) for each k given N=2
         ER2_N3(k)=V0bar*alpha(k)^N(2) + V_adjbar*N(2)*(alpha(k)^(N(2)-1))*(1-alpha(k)); % estimate of the first two terms in eq (2.9) for each k given N=3

     end
    gammasb = reshape(gkb,p*m,1);
    gammasboot(:,b)=gammasb-gammas; % bootstraped empirical distribution of gammas
    
        % Computation of expected payoff when screening level is in [alpha(1),alpha(m-1)]
    for k=1:m-1  
        intmed_N2=trapz(alpha(k:m),ER1_N2(k:m),2); % numerical integration of ER1_N2 from alpha(k) to alpha(m-1) using a trapezoidal rule
        intmed_N3=trapz(alpha(k:m),ER1_N3(k:m),2); % numerical integration of ER1_N3 from alpha(k) to alpha(m-1) using a trapezoidal rule
        ER_N2(:,k)=ER2_N2(k)+intmed_N2; % estimate of the expected payoff given N=2
        ER_N3(:,k)=ER2_N3(k)+intmed_N3; % estimate of the expected payoff given N=3    
    end

    % Expected payoff all alpha in [alpha(1),alpha(m)]
    ER_N2all=[ER_N2,ER2_N2(m)]; % estimate given N=2
    ER_N3all=[ER_N3,ER2_N3(m)]; % estimate given N=3

    % Maximization of expected payoff N=2
    ER_N2max=max(ER_N2all,[],2); % max of expected payoff
    [~,indx]=ismember(ER_N2max,ER_N2all);
    qindx_N2=indx; % index in Alpha that maximize estimated expected payoff
    q_rstar_N2b=alpha(qindx_N2); % optimal screening level
    q_rstar_N2boot(:,b)=q_rstar_N2b-q_rstar_N2; % bootstraped empirical distribution of screening levels
    psi_q_rstar_N2_N2=N(1)*(q_rstar_N2b^(N(1)-1))-(N(1)-1)*(q_rstar_N2b^N(1)); % psi function used in the computation of the optimal reservation price
    psi_q_rstar_N2_N3=N(2)*(q_rstar_N2b^(N(2)-1))-(N(2)-1)*(q_rstar_N2b^N(2)); % psi function used in the computation of the optimal reservation price
    psi_q_rstar_N2 = [ones(L2,1)*psi_q_rstar_N2_N2;ones(L3,1)*psi_q_rstar_N2_N3];
    [~,~,g_q_rstar_N2]=MM(psi_q_rstar_N2, p, wb, x, toler, gammai); % gamma_hat corresponding to the optimal screening level
    
    % Optimal Res Price
    V_rstar_N2b=xbar*g_q_rstar_N2; % optimal reservation price
    V_rstar_N2boot(:,b)=V_rstar_N2b-V_rstar_N2; % bootstraped empirical distribution of optimal reservation price

    % Maximization of expected payoff N=3
    ER_N3max=max(ER_N3all,[],2); % max of expected payoff
    [~,indx]=ismember(ER_N3max,ER_N3all);
    qindx_N3=indx; % index in Alpha that maximize estimated expected payoff
    q_rstar_N3b=alpha(qindx_N3); % optimal screening level
    q_rstar_N3boot(:,b)=q_rstar_N3b-q_rstar_N3; % bootstraped empirical distribution of screening levels
    psi_q_rstar_N3_N2=N(1)*(q_rstar_N3b^(N(1)-1))-(N(1)-1)*(q_rstar_N3b^N(1)); % psi function used in the computation of the optimal reservation price
    psi_q_rstar_N3_N3=N(2)*(q_rstar_N3b^(N(2)-1))-(N(2)-1)*(q_rstar_N3b^N(2)); % psi function used in the computation of the optimal reservation price
    psi_q_rstar_N3 = [ones(L2,1)*psi_q_rstar_N3_N2;ones(L3,1)*psi_q_rstar_N3_N3];
    [~,~,g_q_rstar_N3]=MM(psi_q_rstar_N3, p, wb, x, toler, gammai); % gamma_hat corresponding to the optimal screening level
    
    % Optimal Res Price
    V_rstar_N3b=xbar*g_q_rstar_N3; % optimal reservation price
    V_rstar_N3boot(:,b)=V_rstar_N3b-V_rstar_N3; % bootstraped empirical distribution of optimal reservation price
    
    % Computation of the Seller Expected Revenue under optimal screening level
    ER2_N2_rstar=V0bar*q_rstar_N2b^N(1) + (V_rstar_N2b*N(1))*(q_rstar_N2b^(N(1)-1))*(1-q_rstar_N2b);
    ER2_N3_rstar=V0bar*q_rstar_N3b^N(2) + (V_rstar_N3b*N(2))*(q_rstar_N3b^(N(2)-1))*(1-q_rstar_N3b);
    q2_r=linspace(q_rstar_N2b,0.9,m);
    q3_r=linspace(q_rstar_N3b,0.9,m);
    psiN2_q2=N(1)*(q2_r.^(N(1)-1))-(N(1)-1)*(q2_r.^N(1));
    psiN3_q2=N(2)*(q2_r.^(N(2)-1))-(N(2)-1)*(q2_r.^N(2));
    psiN2_q3=N(1)*(q3_r.^(N(1)-1))-(N(1)-1)*(q3_r.^N(1));
    psiN3_q3=N(2)*(q3_r.^(N(2)-1))-(N(2)-1)*(q3_r.^N(2));

    for k=1:m
        [~,~,g2_r]=MM([ones(L2,1)*psiN2_q2(k);ones(L3,1)*psiN3_q2(k)], p, wb, x, toler, gammai);
        ER1_N2_rstar(:,k)=xbar*g2_r*(N(1)*(N(1)-1)*(q2_r(k)^(N(1)-2))*(1-q2_r(k)));
        [~,~,g3_r]=MM([ones(L2,1)*psiN2_q3(k);ones(L3,1)*psiN3_q3(k)], p, wb, x, toler, gammai);
        ER1_N3_rstar(:,k)=xbar*g3_r*(N(2)*(N(2)-1)*(q3_r(k)^(N(2)-2))*(1-q3_r(k)));
    end
    int2_rstar=trapz(q2_r,ER1_N2_rstar,2); % numerical integration of ER1_N2_rstar using a trapezoidal rule
    int3_rstar=trapz(q3_r,ER1_N3_rstar,2); % numerical integration of ER1_N3_rstar using a trapezoidal rule

    % Expected Revenue under Optimal Reservation Price Policy 
    ER_N2_rstarboot(:,b)=ER2_N2_rstar+int2_rstar - ER_N2_rstar; % bootstraped empirical distribution of the optimal expected payoff
    ER_N3_rstarboot(:,b)=ER2_N3_rstar+int3_rstar - ER_N3_rstar; % bootstraped empirical distribution of the optimal expected payoff

    b
end

% Critical Values
gammas_c025=quantile(gammasboot',0.025); 
gammas_c975=quantile(gammasboot',0.975);
opt_scrN2_c025=quantile(q_rstar_N2boot',0.025);
opt_scrN2_c975=quantile(q_rstar_N2boot',0.975);
opt_scrN3_c025=quantile(q_rstar_N3boot',0.025);
opt_scrN3_c975=quantile(q_rstar_N3boot',0.975);
opt_resN2_c025=quantile(V_rstar_N2boot',0.025);
opt_resN2_c975=quantile(V_rstar_N2boot',0.975);
opt_resN3_c025=quantile(V_rstar_N3boot',0.025);
opt_resN3_c975=quantile(V_rstar_N3boot',0.975);
opt_ERN2_c025=quantile(ER_N2_rstarboot',0.025);
opt_ERN2_c975=quantile(ER_N2_rstarboot',0.975);
opt_ERN3_c025=quantile(ER_N3_rstarboot',0.025);
opt_ERN3_c975=quantile(ER_N3_rstarboot',0.975);


% 95% Confidence Intervals
gammas_CIup=gammas+gammas_c975';
gammas_CIlow=gammas+gammas_c025';
gammas_apv=gk(2,:);
gammas_vol=gk(3,:);
opt_scrN2_CIup=q_rstar_N2+opt_scrN2_c975';
opt_scrN2_CIlow=q_rstar_N2+opt_scrN2_c025';
opt_scrN3_CIup=q_rstar_N3+opt_scrN3_c975';
opt_scrN3_CIlow=q_rstar_N3+opt_scrN3_c025';
opt_resN2_CIup=V_rstar_N2+opt_resN2_c975';
opt_resN2_CIlow=V_rstar_N2+opt_resN2_c025';
opt_resN3_CIup=V_rstar_N3+opt_resN3_c975';
opt_resN3_CIlow=V_rstar_N3+opt_resN3_c025';
opt_ERN2_CIup=ER_N2_rstar+opt_ERN2_c975';
opt_ERN2_CIlow=ER_N2_rstar+opt_ERN2_c025';
opt_ERN3_CIup=ER_N3_rstar+opt_ERN3_c975';
opt_ERN3_CIlow=ER_N3_rstar+opt_ERN3_c025';

% 95% Confidence Intervals for each quantile level
gammas_CIup_k=reshape(gammas_CIup,p,m);
gammas_CIlow_k=reshape(gammas_CIlow,p,m);
CIup_apv=gammas_CIup_k(2,:);
CIlow_apv=gammas_CIlow_k(2,:);
CIup_vol=gammas_CIup_k(3,:);
CIlow_vol=gammas_CIlow_k(3,:);
